#include <stdio.h>
#include <math.h>
#include <cuda.h>
#include <cuda_runtime_api.h>
#define CHECK_ERROR(call) do {       \
    if( cudaSuccess != call) {       \
   fprintf(stderr, "Cuda error in file '%s' in line %i .\n",              \
         __FILE__, __LINE__  );                                           \
         exit(0);                                                         \
    } } while (0)

#define BLOCK_SIZE        512
////////////////////////////////////////////////////////////////////////////////
//
// constantes usadas apenas pelo kernel
//
////////////////////////////////////////////////////////////////////////////////
__constant__ float ALPHA, BETA;

////////////////////////////////////////////////////////////////////////////////
//
// Kernel de 2a no espaco e 2a no tempo
//
////////////////////////////////////////////////////////////////////////////////
__global__ void kernel2H2T(float* C, float* A, float* Aold, unsigned int
pulseOffset, float excitation) {
   unsigned int bx = blockIdx.x;
   unsigned int by = blockIdx.y;
   unsigned int tx = threadIdx.x;
   unsigned int ty = threadIdx.y;
   
   unsigned int width = gridDim.x * blockDim.x;
   unsigned int height = gridDim.y * blockDim.y;
   unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
   unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
   unsigned int kn = y * width + x;
   
   unsigned int k  = ty * blockDim.x + tx;

	extern __shared__ float AS[];

	float tempResult = 0.0f,
	      left	     = 0.0f,
	      right	     = 0.0f,
	      up	        = 0.0f,
	      down	     = 0.0f,
	      old	     = 0.0f,
        
         factor     = 1.0f,
         oldRight1  = 0.0f,
         oldRight2  = 0.0f,
         oldLeft1   = 0.0f,
         oldLeft2   = 0.0f;

	 
	
	AS[k] = A[kn];
  
	__syncthreads();

 
	old = Aold[kn];

	if (tx==0) {							 
		if (bx==0)
      {
			left = 0.0f;
         oldRight1 = Aold[kn+1];
         oldRight2 = Aold[kn+2];
      }
		else
			left = A[kn-1];					 
		right = AS[k+1];				 
      
	} else if (tx==(blockDim.x-1)) {		  
		left = AS[k-1];						   
		if (bx==(gridDim.x-1))
      {
			right = 0.0f;					  
         oldLeft1 =Aold[width*y+x-1];
         oldLeft2 =Aold[width*y+x-2];
      }
		else
			right = A[kn+1];				  

	} else {								   
		left = AS[k-1];
		right = AS[k+1];
	}

	if (ty==0) {						 
		if (by==0)							 
			up = 0.0f;
		else
			up = A[kn-width];    
		down = AS[k+blockDim.x];
	} else if (ty==(blockDim.y-1)) {		   
		up = AS[k-blockDim.x];
		if (by==(gridDim.y-1))				    
			down = 0.0f;
		else
			down = A[kn+width];  
	} else {								   
		up = AS[k-blockDim.x];
		down = AS[k+blockDim.x];
	}

	
	__syncthreads(); 

 
   tempResult = 2.0f * AS[k] - old + ALPHA * (left + right + up + down - 4.0f * AS[k]);
   /*
   if (kn==pulseOffset) 
      tempResult += excitation;
    */ 
 
   //Tem duas bordas absorventes! esquerda e direita
 
   if (tx==0)
   {	
						  
		if (bx==0)
      {
       
         tempResult = AS[k] + right - oldRight1 + factor*(right-AS[k]-(oldRight2-oldRight1));
      }
			
   }else if (tx==(blockDim.x-1))
   {
		
		if (bx==(gridDim.x-1))
      {
         tempResult = AS[k] + left - oldLeft1 - factor*(AS[k]-left-(oldLeft1-oldLeft2));
      }
   }

	C[kn] = tempResult;
}



////////////////////////////////////////////////////////////////////////////////
//
// Variáveis globais da GPU/CPU
//
////////////////////////////////////////////////////////////////////////////////
dim3  dGrid,
      dTthreads;


unsigned int uMemSize = 0,
             uSharedMemSize = 0;

int          iPulseOffset = 0;

float        *fGPUTime  = NULL,  //Time
             *fGPUTime1 = NULL,  //Time - 1
             *fGPUTime2 = NULL,
             *fPtr      = NULL;
             
float *CPUTime;




////////////////////////////////////////////////////////////////////////////////
//
// Função para alocação de recursos na GPU / roda na CPU
//
////////////////////////////////////////////////////////////////////////////////
extern "C" void initGPU(int gpu,
                        unsigned int resX,
                        unsigned int resY,
                        float alpha,
                        float beta)
{

   // Calculando variáveis do CUDA
   unsigned int widthInBlocks,
                heightInBlocks,
                pulsePosX,
                pulsePosY;


  
   CHECK_ERROR(cudaSetDevice(gpu));
   widthInBlocks = (resX / BLOCK_SIZE);
   heightInBlocks = (resY / BLOCK_SIZE);
   //grid(widthInBlocks, heightInBlocks);
   dGrid.x = widthInBlocks;
   dGrid.y = heightInBlocks;
   //threads(BLOCK_SIZE, BLOCK_SIZE);
   dTthreads.x = BLOCK_SIZE;
   dTthreads.y = BLOCK_SIZE;

   // Calculando offset do pulso

   pulsePosX   = (unsigned int)floor(0.5f * resX);
   pulsePosY   = (unsigned int)floor(0.0f * resY);
   iPulseOffset =  pulsePosY * resX + pulsePosX;

   // Alocação
   uSharedMemSize = dTthreads.x * dTthreads.y * sizeof(float);
   uMemSize = resX * resY * sizeof(float);
   fprintf(stdout, "\nMemória usada na GPU [%d]bytes", uMemSize*3);

   CPUTime = (float *) calloc(uMemSize, sizeof(float));
    
   
   int x, y, p;
   float X,  Y, ret;
   float W = (float) (dGrid.x * dTthreads.x);
   float H = (float) (dGrid.y * dTthreads.y);
   
   //Iniciakuza pulso
   for (x = 0; x < (dGrid.x * dTthreads.x); x++)
   {
      for (y = 0; y < (dGrid.y * dTthreads.y); y++)
      {
         X = (x - (0.5f * W)) * (x - (0.5f * H));
         Y = (y - (0.5f * W)) * (y - (0.5f * H));
         ret = exp ( -(X + Y) / (0.0000250f  * H  * W)  ); //0.000075
         p = y * (dGrid.x * dTthreads.x) + x;
         CPUTime[p] = ret;
      }
   }



   CHECK_ERROR(cudaMalloc((void**) &fGPUTime, uMemSize));
   CHECK_ERROR(cudaMalloc((void**) &fGPUTime1, uMemSize));
   CHECK_ERROR(cudaMalloc((void**) &fGPUTime2,  uMemSize));

   CHECK_ERROR(cudaMemcpy(fGPUTime1, CPUTime, uMemSize, cudaMemcpyHostToDevice));
   //CHECK_ERROR(cudaMemcpy(fGPUTime2, CPUTime, uMemSize, cudaMemcpyHostToDevice));
   free(CPUTime);
   
   
   //CHECK_ERROR(cudaMemset(fGPUTime1,  0, uMemSize));
   CHECK_ERROR(cudaMemset(fGPUTime2,  0, uMemSize));

   //copia das constantes
   CHECK_ERROR(cudaMemcpyToSymbol("ALPHA", &alpha, sizeof(float)));
   CHECK_ERROR(cudaMemcpyToSymbol("BETA", &beta, sizeof(float)));
//   CHECK_ERROR(cudaMemcpyToSymbol("DELTAT", &deltaT, sizeof(float)));
      


   //Aloca o vetor na CPU se for salvar em arquivo
  
   fprintf(stdout, "\nDiscretização: %d x %d", resX, resY);
   fprintf(stdout, "\n         Grid:  %d x %d", widthInBlocks, heightInBlocks);
   fprintf(stdout, "\n        Block: %d x %d", BLOCK_SIZE, BLOCK_SIZE);
   fprintf(stdout, "\n");
   fflush(stdout);
}

////////////////////////////////////////////////////////////////////////////////
//
// Função para chamar kernel 
//
////////////////////////////////////////////////////////////////////////////////
extern "C" void runWave(float excitation, float *d)
{
  
   
  
   kernel2H2T<<<dGrid, dTthreads, uSharedMemSize>>> (fGPUTime, fGPUTime1, fGPUTime2, iPulseOffset, excitation);

   CHECK_ERROR(cudaThreadSynchronize()); //sincroniza threads
   

   CHECK_ERROR(cudaMemcpy(d, fGPUTime, uMemSize, cudaMemcpyDeviceToHost));
   
   fGPUTime2 = fGPUTime1;
   fGPUTime1 = fGPUTime;
   fGPUTime = fGPUTime2;
 
}


////////////////////////////////////////////////////////////////////////////////
//
// Função para alocação de recursos na GPU / roda na CPU
//
////////////////////////////////////////////////////////////////////////////////
extern "C" void freeGPU(void)
{
   cudaFree(fGPUTime);
   cudaFree(fGPUTime1);
   cudaFree(fGPUTime2);
}

